
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# DISCLAIMER AND GENERAL INFORMATION
#
# File: 01_main_analysis.R
# Purpose: Main data analysis
# Date: June 2025
# Data: pulled through 00_data_prep.R
#
# See 00_data_prep.R for technical disclaimer on R versions
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Load data
source("00_data_prep.R")


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (A) Descriptive results ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


# General population means for all four levels
mean(full$dv_climate_intl_dummy[full$survey=="Prolific"])
mean(full$dv_climate_nat_dummy[full$survey=="Prolific"])
mean(full$dv_climate_subnat_dummy[full$survey=="Prolific"])
mean(full$dv_climate_local_dummy[full$survey=="Prolific"])

# Descriptive information on mitigation/adaptation preferences

table(full$sample, full$mitig_v_adapt)

N <- table(full$sample) - table(full$mitig_v_adapt,full$sample, useNA="ifany")[4,]
df <- data.frame(matrix(rep(N,3),nrow=3, ncol=4, byrow = T))
m <- table(full$mitig_v_adapt,full$sample)/df
colnames(m) <- names(N)


table(full$sample, full$citizen)
N <- table(full$sample) - table(full$sample, full$citizen, useNA="ifany")[,5]
df <- data.frame(matrix(rep(N,4),nrow=4, ncol=4))
m <- table(full$sample, full$citizen)[,-c(5)]/df
rownames(m) <- names(N)
colnames(m) <- colnames(table(full$sample, full$citizen)[,-c(5)])
m


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (B) Main results: Multinomial models with controls ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Define control variables
ctrls <- paste(c("factor(female)", "factor(middleage)", "factor(old)", "factor(educ_cat)", "factor(pol_left)", "factor(pol_right)"), collapse="+")

# Estimate model
m2 <- multinom(paste("factor(effect_cc_1st) ~ sample", ctrls, sep="+"), data=full)

# Get predicted probabilities
d2 <- plot_predictions(m2, type="probs", by="sample", draw=FALSE) %>%
  mutate(group=as.numeric(group)) %>% 
  arrange(group,sample) %>%
  mutate(rowid=rep(seq(1,4,1),4))

# Create plotting data
dt <- avg_predictions(m2, type="probs", by="sample", hypothesis = ~ pairwise) %>%
filter(hypothesis=="(1 Scotland) - (1 Northern England)" |
       hypothesis=="(1 Wales) - (1 Northern England)" |   
       hypothesis=="(2 Scotland) - (2 Northern England)" |
       hypothesis=="(2 Wales) - (2 Northern England)" |
       hypothesis=="(3 Scotland) - (3 Northern England)" |
       hypothesis=="(3 Wales) - (3 Northern England)" |
       hypothesis=="(4 Scotland) - (4 Northern England)" |
       hypothesis=="(4 Wales) - (4 Northern England)") %>%
mutate(rowid=case_when(grepl("Scotland",hypothesis) ~ 3,
    grepl("Wales",hypothesis) ~ 4),
    group=case_when(grepl("1",hypothesis) ~ 1,
    grepl("2",hypothesis) ~ 2,
    grepl("3",hypothesis) ~ 3,
    grepl("4",hypothesis) ~ 4),
    p5=case_when(p.value<0.05 ~ "yes", .default = "no"),
    p10=case_when(p.value<0.1 ~ "yes", .default = "no")) %>%
select(rowid,group,hypothesis,p5,p10) %>%
right_join(d2,by = join_by("rowid","group")) %>%
arrange(group, rowid)

# Change facet labels
group.labs = c("1"="Local level", "2"="Devolved nation level", "3"="National level", "4"="International level")

# Change group labels 
dt <- dt %>%
  mutate(sample=case_when(sample=="Northern England" ~ "Yorkshire and Cumbria",
                          sample!="Northern England" ~ sample,
                          .default=NA)) %>%
  mutate(sample=factor(sample, levels=c("General population", "Yorkshire and Cumbria", "Scotland", "Wales"), 
                       ordered=FALSE)) %>%
  as_tibble()


# Create results plot
p <- ggplot(data=dt,
            aes(x=estimate, y=sample, group=group, 
                         color=p5)) +
  geom_rect(data=NULL,aes(xmin=0,xmax=Inf,ymin=.5,ymax=1.5),
            fill=cols[1], alpha=.1, color=NA) +
  geom_point() + 
  geom_vline(xintercept=0, linetype="dashed") +
  geom_linerangeh(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.8)) +
  geom_text(aes(label = ifelse(is.na(p5)==FALSE & p5=="yes", "*", ""), group=group),
            position = position_dodge(width = .8), vjust=.1) + 
  facet_wrap(~group, labeller = labeller(group = group.labs)) +
  ylab("") +
  xlab("Predicted probability for governance level") + 
  scale_color_manual(values=c("black",cols[2]), na.value="black") + 
  theme_classic() +
  theme(legend.position="none")
p

ggsave(p, file="./plots/main_results.pdf", width=6, height=4)

# # Produce main output file in jpg
# jpeg(file="./plots/main_results.jpg", width=6, height=4, units="in", res=600)
# p
# dev.off()

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#                         END OF FILE
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


